************************************************************
* Predictive Modeling of Women's Mental Health:
* Financial Confidence as a Humanitarian Determinant in Bangladesh
*
* Conference:
* 2025 IEEE International Women in Engineering (WIE) Conference
* on Electrical and Computer Engineering (WIECON ECE)
* Cox's Bazar, Bangladesh
*
* Authors:
* Asma Binte Afzal
* Department of Environmental and Occupational Health
* Indiana University Bloomington, Indiana, USA
* ORCID 0000 0003 4318 376X
*
* Muhammad Hassan Bin Afzal
* Department of Applied Statistics and Data Science
* Jahangirnagar University, Dhaka, Bangladesh
* ORCID 0000 0001 8192 0885
*
* Core research question
* To what extent does women's perceived financial confidence
* in accessing healthcare predict anxiety and depression, and
* how are these associations shaped by education, digital
* inclusion, and community context in Bangladesh
*
* This script assumes
* 1 The custom analytic dataset is already loaded in memory
* 2 The dataset is fully anonymized and includes only derived
*   variables needed for this paper
*   gad7_ord4 phq9_ord5 finconf_med autonomy_say bdlift
*   education wealth age is_working residence bank_account
*   wt cluster_id_anon analytic_sample
************************************************************

************************************************************
* 1 Analytic sample restriction
************************************************************

* Keep analytic sample if the flag exists
capture confirm variable analytic_sample
if _rc==0 {
    keep if analytic_sample==1
}

************************************************************
* 2 Basic labels and quick checks
************************************************************

* Make sure yes or no labels exist
capture label define yesno 0 "No" 1 "Yes"
capture label values is_working yesno
capture label values bank_account yesno

* Tabulate core variables
tab gad7_ord4
tab phq9_ord5
tab finconf_med
tab bdlift
tab education
tab wealth
tab age
tab is_working
tab residence
tab bank_account

summ gad7_ord4 phq9_ord5 wt

* Simple histograms for outcome distributions
histogram gad7_ord4, discrete freq ///
    title("Distribution of GAD 7 categories")
histogram phq9_ord5, discrete freq ///
    title("Distribution of PHQ 9 categories")

************************************************************
* 3 Survey design declaration for descriptive checks
*   Multilevel models below use cluster random effects with weights
************************************************************

svyset cluster_id_anon [pweight = wt]

************************************************************
* 4 Descriptive statistics with survey design
************************************************************

* Outcome distributions
svy: tab gad7_ord4, percent
svy: tab phq9_ord5, percent

* Key predictors and controls
svy: tab finconf_med, percent
svy: tab bdlift, percent
svy: tab education, percent
svy: tab wealth, percent
svy: tab age, percent
svy: tab is_working, percent
svy: tab residence, percent
svy: tab bank_account, percent

************************************************************
* 5 Correlation and multicollinearity diagnostics
************************************************************

* Pairwise correlations among predictors
pwcorr finconf_med autonomy_say bdlift education wealth age ///
       is_working residence bank_account, sig star(0.05)

* Multicollinearity check with collin
capture which collin
if _rc ssc install collin

collin finconf_med autonomy_say bdlift education wealth age ///
       is_working residence bank_account

************************************************************
* 6 Multilevel ordinal logistic models (main results)
*   Outcomes gad7_ord4 and phq9_ord5
************************************************************

* GAD 7 bivariate model
meologit gad7_ord4 i.finconf_med [pweight = wt] ///
    || cluster_id_anon:, or

estimates store GAD1

* GAD 7 full model
meologit gad7_ord4 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

estimates store GAD2

* PHQ 9 bivariate model
meologit phq9_ord5 i.finconf_med [pweight = wt] ///
    || cluster_id_anon:, or

estimates store PHQ1

* PHQ 9 full model
meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

estimates store PHQ2

* Intraclass correlation for GAD 7 full model
di "Intraclass correlation for GAD 7 full model"
estat icc

* Re run PHQ 9 full model to get ICC for that outcome
meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

di "Intraclass correlation for PHQ 9 full model"
estat icc

* Optional side by side comparison table if esttab is installed
capture which esttab
if !_rc {
    esttab GAD1 GAD2 PHQ1 PHQ2, eform ///
        cells(b(star fmt(3)) se(par fmt(3))) ///
        stats(N ll bic aic, fmt(0 3 3 3) ///
              labels("Obs" "Log Likelihood" "BIC" "AIC")) ///
        legend collabels(none) ///
        starlevels(* 0.05 ** 0.01 *** 0.001) ///
        mtitles("GAD7 Bivariate" "GAD7 Full" ///
                "PHQ9 Bivariate" "PHQ9 Full")
}

************************************************************
* 7 Calibration plots for full multilevel models
*   GAD 7 and PHQ 9
************************************************************

* GAD 7 full model for calibration
meologit gad7_ord4 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

predict phat_gad if e(sample), pr

* Create deciles of predicted probability
xtile decile_gad = phat_gad, n(10)

* Observed and predicted by decile
egen mean_phat_gad     = mean(phat_gad),    by(decile_gad)
egen mean_observed_gad = mean(gad7_ord4),   by(decile_gad)

twoway (scatter mean_observed_gad mean_phat_gad, msymbol(circle)) ///
       (line mean_observed_gad mean_phat_gad, sort), ///
       title("Calibration GAD 7 full model") ///
       xtitle("Predicted probability") ///
       ytitle("Observed mean category") ///
       legend(off)

graph save cal_gad7.gph, replace

* PHQ 9 full model for calibration
meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

predict phat_phq if e(sample), pr

xtile decile_phq = phat_phq, n(10)

egen mean_phat_phq     = mean(phat_phq),    by(decile_phq)
egen mean_observed_phq = mean(phq9_ord5),   by(decile_phq)

twoway (scatter mean_observed_phq mean_phat_phq, msymbol(circle)) ///
       (line mean_observed_phq mean_phat_phq, sort), ///
       title("Calibration PHQ 9 full model") ///
       xtitle("Predicted probability") ///
       ytitle("Observed mean category") ///
       legend(off)

graph save cal_phq9.gph, replace

graph combine cal_gad7.gph cal_phq9.gph, col(2) ///
    title("Calibration plots financial confidence and mental health")

************************************************************
* 8 Predictive margins and interaction plots
*   Focus mild category outcome(1)
************************************************************

* Education × financial confidence mild anxiety GAD 7
meologit gad7_ord4 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins education, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Education and financial confidence Mild anxiety GAD 7") ///
    ytitle("Pr Mild anxiety") ///
    xtitle("Education level") ///
    by(finconf_med) ///
    name(gad_mild_finconf_edu, replace)

* Education × financial confidence mild depression PHQ 9
meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins education, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Education and financial confidence Mild depression PHQ 9") ///
    ytitle("Pr Mild depression") ///
    xtitle("Education level") ///
    by(finconf_med) ///
    name(phq_mild_finconf_edu, replace)

graph combine gad_mild_finconf_edu phq_mild_finconf_edu, col(2) ///
    title("Education and financial confidence Mild mental health")

* Age × financial confidence
meologit gad7_ord4 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins age, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Age and financial confidence Mild anxiety GAD 7") ///
    ytitle("Pr Mild anxiety") ///
    xtitle("Age groups") ///
    by(finconf_med) ///
    name(gad_mild_finconf_age, replace)

meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins age, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Age and financial confidence Mild depression PHQ 9") ///
    ytitle("Pr Mild depression") ///
    xtitle("Age groups") ///
    by(finconf_med) ///
    name(phq_mild_finconf_age, replace)

graph combine gad_mild_finconf_age phq_mild_finconf_age, col(2) ///
    title("Age and financial confidence Mild mental health")

* Wealth × financial confidence
meologit gad7_ord4 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins wealth, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Wealth and financial confidence Mild anxiety GAD 7") ///
    ytitle("Pr Mild anxiety") ///
    xtitle("Wealth quintile") ///
    by(finconf_med) ///
    name(gad_mild_finconf_wealth, replace)

meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins wealth, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Wealth and financial confidence Mild depression PHQ 9") ///
    ytitle("Pr Mild depression") ///
    xtitle("Wealth quintile") ///
    by(finconf_med) ///
    name(phq_mild_finconf_wealth, replace)

graph combine gad_mild_finconf_wealth phq_mild_finconf_wealth, col(2) ///
    title("Wealth and financial confidence Mild mental health")

* Employment × financial confidence
meologit gad7_ord4 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins is_working, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Employment and financial confidence Mild anxiety GAD 7") ///
    ytitle("Pr Mild anxiety") ///
    xtitle("Employment status") ///
    xlabel(0 "Not working" 1 "Working") ///
    by(finconf_med) ///
    name(gad_mild_finconf_work, replace)

meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins is_working, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Employment and financial confidence Mild depression PHQ 9") ///
    ytitle("Pr Mild depression") ///
    xtitle("Employment status") ///
    xlabel(0 "Not working" 1 "Working") ///
    by(finconf_med) ///
    name(phq_mild_finconf_work, replace)

graph combine gad_mild_finconf_work phq_mild_finconf_work, col(2) ///
    title("Employment and financial confidence Mild mental health")

* Bank account × financial confidence
meologit gad7_ord4 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins bank_account, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Bank account and financial confidence Mild anxiety GAD 7") ///
    ytitle("Pr Mild anxiety") ///
    xtitle("Bank account") ///
    xlabel(0 "No account" 1 "Has account") ///
    by(finconf_med) ///
    name(gad_mild_finconf_bank, replace)

meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins bank_account, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("Bank account and financial confidence Mild depression PHQ 9") ///
    ytitle("Pr Mild depression") ///
    xtitle("Bank account") ///
    xlabel(0 "No account" 1 "Has account") ///
    by(finconf_med) ///
    name(phq_mild_finconf_bank, replace)

graph combine gad_mild_finconf_bank phq_mild_finconf_bank, col(2) ///
    title("Bank account and financial confidence Mild mental health")

* BD LIFT × financial confidence
meologit gad7_ord4 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins bdlift, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("BD LIFT and financial confidence Mild anxiety GAD 7") ///
    ytitle("Pr Mild anxiety") ///
    xtitle("BD LIFT index") ///
    by(finconf_med) ///
    name(gad_mild_finconf_bdlift, replace)

meologit phq9_ord5 i.finconf_med i.autonomy_say i.bdlift ///
    i.education i.residence i.bank_account i.wealth i.age ///
    i.is_working [pweight = wt] ///
    || cluster_id_anon:, or

margins bdlift, at(finconf_med=(0 1)) predict(outcome(1))
marginsplot, ///
    title("BD LIFT and financial confidence Mild depression PHQ 9") ///
    ytitle("Pr Mild depression") ///
    xtitle("BD LIFT index") ///
    by(finconf_med) ///
    name(phq_mild_finconf_bdlift, replace)

graph combine gad_mild_finconf_bdlift phq_mild_finconf_bdlift, col(2) ///
    title("BD LIFT and financial confidence Mild mental health")

************************************************************
* 9 Proportional odds checks Brant and gologit2
*   These are single level checks for the parallel lines assumption
************************************************************

* Install brant and gologit2 if needed
capture which brant
if _rc ssc install brant

capture which gologit2
if _rc ssc install gologit2

* GAD 7 ordered logit with Brant test
ologit gad7_ord4 finconf_med autonomy_say bdlift ///
       education wealth age is_working

brant, detail

* PHQ 9 ordered logit with Brant test
ologit phq9_ord5 finconf_med autonomy_say bdlift ///
       education wealth age is_working

brant, detail

* Generalized ordered logit for GAD 7
gologit2 gad7_ord4 finconf_med autonomy_say bdlift ///
         education wealth age is_working, autofit

* Generalized ordered logit for PHQ 9
gologit2 phq9_ord5 finconf_med autonomy_say bdlift ///
         education wealth age is_working, autofit

************************************************************
* End of replication script
************************************************************

